#include <iostream>
#include <cstdlib>
#include <vector>
#include <omp.h>
#include <sys/resource.h>
#include <ctime>
#include <cmath>
#include <algorithm>

using namespace std;


void printElapsed (char *desc, struct timeval *start, struct timeval *end, int niters )
{
	struct timeval elapsed;
	if (start -> tv_usec > end -> tv_usec)
	{
		end -> tv_usec += 1000000;
		end -> tv_sec--;
	}

	elapsed.tv_usec = end -> tv_usec - start -> tv_usec;
	elapsed.tv_sec = end -> tv_sec - start -> tv_sec;

	printf("%s\n total elapsed time = %1d (usec)\n",
			 desc, (elapsed.tv_sec*1000000 + elapsed.tv_usec)/niters);

}/* printElapsed() */

double realrand()
{
	return rand () / (double) (RAND_MAX);
}

/*The computing domain is bounded by [0, 1]*/
void generate_sites (unsigned * tags, double * sites)
{
#pragma omp parallel shared (tags, sites)
	{
		int thread_id = omp_get_thread_num();
		for (int i = tags[thread_id] ; i < tags[thread_id + 1]; ++i)
		{
			sites[i] = realrand();
		}
	}
}

void training(unsigned * tags, double * sites, 
			  unsigned * weights, unsigned * local_min, 
			  unsigned  num_threads, unsigned num_iterations)
{
	double seeds = 0.0;
	for (int loops = 0; loops < num_iterations; ++loops)
	{
		seeds = realrand();
#pragma omp parallel shared (tags, sites, weights, seeds)
{
	    // Find minimum on this processor
	    unsigned thread_id = omp_get_thread_num();
		unsigned pos = tags[thread_id] ;
		for (int i = tags[thread_id] + 1; i < tags[thread_id + 1] ; ++i)
		{
			if (abs(sites[i] - seeds ) < abs(sites[pos] - seeds))
			{
				pos = i; 
			}
		}
		local_min [thread_id] = pos; 
}

        // Balanced Tree technique to extract the global minimum


    unsigned length = num_threads;
    unsigned * local_min_bk = new unsigned [16];
    while (length > 1)
    {
//	    for (int i = 0; i < length; ++i)
//	    {
//		    cout<<local_min[i]<<' ';
//	    }
//	    cout<<endl;

#pragma omp parallel shared (seeds, sites, local_min, local_min_bk )
{
	    unsigned thread_id = omp_get_thread_num();
	
	    if (thread_id % 2 == 0 && thread_id < length )
	    {            
			if (thread_id + 1 < length)
			{
				if (abs (sites[local_min[thread_id]] - seeds) < abs(sites [local_min[thread_id+1]] - seeds))
				{
					local_min_bk[thread_id/2] = local_min[thread_id] ;
				}
				else
				{
					local_min_bk[thread_id/2] = local_min[thread_id + 1] ;
				}
			}
			else
			{
				local_min_bk[thread_id/2] = local_min[thread_id];
			}
	    }
}
#pragma omp parallel shared(local_min, local_min_bk)
{
	unsigned thread_id = omp_get_thread_num();
		local_min[thread_id] = local_min_bk[thread_id];
}
	    length /= 2;

    }

    unsigned global_min_pos = local_min[0];
    sites[global_min_pos] = (sites[global_min_pos] * weights[global_min_pos]
						 + seeds ) / (weights[global_min_pos] + 1);
    weights[global_min_pos] ++;

/*
        //Find the global minimum 
        //cout<<"SEEDS"<< seeds <<endl;
		unsigned idx = 0;
		for (int i = 1; i < num_threads; ++i)
		{
			//cout<<local_min[i]<<endl;
			if (abs (sites[local_min[i]] - seeds) < abs(sites [local_min[idx]] - seeds))
			{
				idx = i;
			}
		}
		unsigned global_min_pos = local_min[idx];
		sites[global_min_pos] = (sites[global_min_pos] * weights[global_min_pos]
								 + seeds ) / (weights[global_min_pos] + 1);
		weights[global_min_pos] ++;
		//cout<<global_min_pos<<endl;
		*/
    }

}

void divide_group (unsigned size, unsigned num_threads, unsigned * tags)
{
	unsigned group_size ;
	int thread_id, d;
	unsigned start , end ;
#pragma omp parallel private (group_size, thread_id, d, start, end) \
shared(tags, size, num_threads)
{

    group_size = size / num_threads ;

	thread_id = omp_get_thread_num();
	d = size % num_threads;
	start = 0;
	end = 0;
	if (thread_id < d) 
	{
		start = thread_id * (group_size + 1 ); 
	}
	else
	{
		start = d * (group_size + 1)  + (thread_id - d) * (group_size); 
	}
	tags[thread_id] = start;
	//cout<<"START TAG PO"<<start<<endl;
}
}
int main( int argc, char * argv[])
{

	int num_sites = 0;                  /* num of sites */
	int num_threads = 0;               /* num of threads */


	struct timeval tmstart, tmend;
	struct timeval genstart, genend;
	struct timezone tzp;

	/*----------------------------------------------------------------
	 * Check arguments
	 *---------------------------------------------------------------*/
	if (argc < 3)
	{
		printf("Usage : %s [num of sites] [num of iterations]\n", argv[0]);
		return 0;
	}
	num_threads = omp_get_max_threads();
	if (num_threads > 16)
	{
		printf("Usage : %s [numints]\n - threads overflow. Maxium number of threads is 16.", argv[0]);
		return 0;
	}

	num_sites = atoi(argv[1]);
	printf("\nnum_threads = %d, num_sites = %d\nExcuting %s ...\n", omp_get_max_threads(), num_sites, argv[0]);


	//Step 1. Divide group.

	unsigned  * tags = new unsigned [num_threads+1];
	tags[num_threads] = num_sites;
	gettimeofday(&genstart, &tzp);
	divide_group (num_sites, num_threads, tags);
//	for (int i = 0;  i < num_threads+ 1; ++i)
//	{
//		cout<<tags[i]<<' ';
//	}
//	cout<<endl;
    gettimeofday(&genend, &tzp);
    printElapsed("\n Step 1  Divide Group  ", &genstart, &genend, 1);

	
	//Step 2. Choose initial sites ;
	double * sites = new double [num_sites];
	gettimeofday(&genstart, &tzp) ;
	generate_sites (tags, sites) ;

//	for (int i = 0; i < num_sites; ++i)
//	{
//		cout<<i<<' '<<sites[i]<<endl;
//	}
//	cout<<endl;

	gettimeofday(&genend, &tzp);
	printElapsed("\n Step 2 Initial data generation ", &genstart, &genend, 1);


	//Step 3. Training the sites ;
	unsigned * weights = new unsigned [num_sites];
	for (int i = 0; i < num_sites; ++i ) weights [i] = 1;


	unsigned * local_min = new unsigned [num_threads];
	unsigned num_iterations = atoi(argv[2]); 
	gettimeofday(&genstart, &tzp);
	training(tags, sites, weights, local_min, num_threads, num_iterations) ;

    gettimeofday(&genend, &tzp);
    printElapsed("\n Step 3  Training the data  ", &genstart, &genend, 1);

	sort (sites, sites + num_sites);
//	for (int i = 0; i < num_sites; ++i)
//	{
//		cout<<i << ' '<<sites[i]<<' '<<weights[i]<<endl;
//	}
//	cout<<endl;

	//Step 4. Release all the memories.
	delete []sites ;
	delete []weights ;
	delete []tags ;
	delete []local_min ;
	return 0;
}
